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We study the relaxation dynamics of a coarse-grained polymer chain at different degrees of stretch- 
ing by both analytical means and numerical simulations. The macromolecule is modelled as a string 
of beads, connected by anharmonic springs, subject to a tensile force applied at the end monomer of 
the chain while the other end is fixed at the origin of coordinates. The impact of bond non-linearity 
on the relaxation dynamics of the polymer at different degrees of stretching is treated analytically 
within the Gaussian self-consistent approach (GSC) and then compared to simulation results derived 
from two different methods: Monte-Carlo (MC) and Molecular Dynamics (MD). 

At low and medium degrees of chain elongation we find good agreement between GSC predictions 
and the Monte-Carlo simulations. However, for strongly stretched chains the MD method, which 
takes into account inertial effects, reveals two important aspects of the nonlinear interaction between 
qq • monomers: (i) a coupling and energy transfer between the damped, oscillatory normal modes of the 

chain, and (ii) the appearance of non-vanishing contributions of a continuum of frequencies around 
the characteristic modes in the power spectrum of the normal mode correlation functions. 
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I. INTRODUCTION 



The dynamics of a linear polymer was first studied and solved by Rouse in his seminal paper— on the relaxation 
behavior of a phantom (ideal) chain of volumclcss beads connected by harmonic (linear) springs. This model, which 
reveals the transcendental features of a linear chain, gives all characteristic dynamic values for nearly free chains 
whereby the assumed linearity of the model provides a good description of the process under investigation. Since 
this simplified model fails to retain the physics when the chain is subjected to strong tensile forces (large stretching), 
^ ■ several approaches have been suggested which take into account the finite extensibility of the chain. Most of them 
are based on modifying the underlying force between monomers, whose nonlinear nature manifests itself for high 
extensions. Among the earlier studies that have been tried to understand the dynamics of highly stretched polymers 
one should mention the work of Pincus^ which considered internal modes of a strongly stretched chain using the 
, blob model^. Based on this work, Marciano and Brochart-Wyarl4, studied a single chain stretched by a force / 
applied at the free end and fixed with the other end at the origin. They demonstrated that stretched chains can be 
thought as Rouse chains of impenetrable blobs where all the blobs are identical. In the long wave-length limit chain 
dynamics is described by renormalized Rouse modes of mode number p with a dispersion relation for the corresponding 
relaxation times t p cx p~ 2 . More recently, the work of Hatfield and Quake^ presented theoretical calculations and 
y—i . computer simulations to estimate the effects of tension and hydrodynamics on r, the fundamental relaxation time of 
' the polymer. 

Apparently none of these studies, however, seems to have reached a complete and definite description of the problem, 
in spite of its fundamental nature. In view of many experiments and technological applications, the dynamics of 
highly stretched chains results fundamental in studying bond rupture of filled and unfilled elastomers. So far, there 
no theory based on molecular models for the fracture of materials exists. Most of the existing theories are based 
on phenomenological ideas, which consider directly the material parameters^^ ^ 10 ' 11 ! 12 . In this work, the main 
emphasis is drawn on developing molecular pictures for breaking chains, unfilled and filled elastomers. It turns out 
quickly that the basic problems concerning the behavior of strongly stretched chains are not treated in such a way, 
that a systematic study on the chain breaking mechanism can be carried out. Therefore we concentrated ourselves 
basically to study the dynamics of highly stretched chains. To do so, there are two different ways to tackle this 
problem. The most obvious possibility is to introduce constraints, where the lengths of the chain segments are fixed 
at any times during the dynamical processes. Such a theory requires the introduction of Lagrange multipliers. The 
theory can be formulated on an exact way, but nevertheless the solution needs many approximations, which are not 
always easy to control. The corresponding stochastic equations are highly non-linear and do not allow the complete 
determination of the Langrange multipliers^. 

To be more specific, we state the problem which we are going to discuss below in some more details. Polymer 
dynamics is usually studied by the so-called Rouse equation, which has its origin in the corresponding Edwards 
Hamiltonian for Gaussian Chains. In discrete notation Gaussian chains are described by a Gaussian chain statistics, 
which is defined by the Hamiltonian 
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where k B T is the thermal energy, b; = — rj_i the individual bond vectors of a chain defined by a connected set of 
position vectors {r^} and 6§ their mean value. The Rouse equation, which is the Langevin equation for the monomer 
position vector is simply defined by the equation of motion 

dR(s,t) 3k B T d 2 R( S ,t) 
t°— t &g c^ =h(s ' i} (2) 

where has been represented by continuous space curves R(s, i) and s is the arc length variable running from to L 
along chain contour. In last equation, the first term corresponds to the frictional force, £o being the friction coefficient 
and h(s, t) is the stochastic force, which is taken as a white noise. The second term is nothing but the elastic force, 
and is easily derived by the gradient of the Hamiltonian H. The drawback of this simple stochastic equation is, that 
the bond vectors can be stretched to infinity and thus Gaussian chains overstretch. 

This drawback is normally repaired by the use of rigid constraints, which fix the bond length strictly to |bj| = 
|rj — = 60, but this yields to all the problems discussed above. The alternative possibility is the use of potentials, 
which do not allow the chain to extend to infinity at high external forces instead of hard constraints. Although this 
corresponds to a very different approach the problems become to some extend more feasible. As an example we 
mention the FENE - potential^, which is widely used in computer simulations, to prevent the overstretching of the 
chains and avoiding so unphysical conformations. To do so we have proposed a modified chain potential according to 
the FENE - potential, where the chain segments can only stretch to a maximum value, i.e., 

where we have chosen here just for the maximum extension of the bonds the mean bond size from above (in this eq. 
this is done for illustration, for the more specific calculations and simulations we have chosen more general cases). The 
distribution function for the end-to-end vector contains Bessel functions that makes it not very suitable for analytical 
computations. Nevertheless in a crude approximation it can be shown, that its asymptotics is of the form 

G(R, N) cx exp (^ (1 _ W( ^ q))2)1/2 ) (4) 

which shows that the end vector distribution function G(R, N) tends to zero when the chain is close to its maximum 
extension R = Nbg. 

Concerning the dynamics we can also extract some features of the FENE - type potential. For small elongations, 
|bj| << bo the FENE - Hamiltonian yields back the Gaussian chain and so, in the dynamics the classical Rouse 
equation. On the other hand it provides the basis for a non-linear dynamical polymer model, which can be summarized 
in the (approximate) non-linear stochastic equation, 

dt % 1 1 f an(s,t) Y 

1 l? { 9s ) 

Our intention in the present work is to examine the fundamental aspect of bond nonlinearity of the chain in the high 
stretched limit and its impact on the polymer relaxation dynamics with analytical methods and simulation techniques. 
To this end we first employ the Gaussian self-consistent approach. The method has been extensively used (see e.gJ^) 
and appears to provide an adequate approach to study this type of problems. In addition, we run Monte-Carlo (MC) 
simulations to compute and compare the dynamical behavior with that, predicted by the analytical approach. It is 
worth noting that the MC scheme does not take into account mass in the model. In this sense, we want to know 
to what extent this is a good description of the problem. To answer this question, we also performed Molecular 
Dynamics (MD) simulations since at high degrees of chain extension the effects of mass and inertia may not be 
neglected; therefore we have to modify the equation of motion as follows: 



3 2 R(s,t) dK(s,t) 3k B T 
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(6) 



The extension by the inertia term is to first instance not obvious, though it is intuitive. The dynamics of strongly 
stretched chains will always get vibrational modes, whenever non-linear springs are connecting the beads. In contrast 



to the Gaussian chains, which can be stretched to infinity, the springs take at large stretching forces more and more 
energy, which will manifest itself in travelling waves along the longitudinal direction of the stretched string. Of course 
this FENE Rouse model is also only an approximation, since rigid constraints cannot be described rigorously by any 
(soft) potential. Nevertheless it removes the basic non-physical features from the Gaussian chain model. The chain 
possesses a finite extensibility and the distribution function does not contain unphysical (overstretched) states. 
Form this equation we can expect several new features compared to the Rouse equation. 

1. We can expect diverging time scales as long as the chain becomes strongly extended, i.e., close to the maximum 
deformation, where only chain fluctuations close to the maximum stretching of the individual bonds matter. 
These diverging time scales will be imposed by the singularity in the stretching force term. 

2. We can expect a "transition" from relaxational dynamics to a travelling wave dynamics. At low deformation 
states the dynamics is mainly ruled by fluctuations where the FENE nature does not play a crucial role. Thus 
at low deformations we readily expect a Rouse type dynamics. At larger stretching ratios, we have basically a 
linear chain (like in one dimensional mono atomic solids) and "acoustic phonon modes" can be expected, if the 
springs were harmonic. 

3. Since the springs are highly anharmonic we can expect non-trivial features for the chain dynamics. At high 
stretching ratios, i.e., at large stretching forces, we might have strong effects from the non-linearity of the poten- 
tial. At small deformations the Rouse equation is still a valid and useful approximation, at large deformations, 
however, the "Rouse modes" will interact with each other, they will exchange energy and damped solitary wave 
excitations will appear. 

However for the main conclusions the form and power of the singularity in the modified force term in eq.([l]) is not 
very important. The main issue is that a singularity at full extension of the bonds, i.e, the chain appears. Another 
aspect is that we have so far ignored the orientation of the chain upon stretching, which will simplify the problem 
slightly since we will split the isotropic variable R(s) into its direction perpendicular and parallel to the force. It is 
then physically obvious that both components will respond differently to the strong stretching force. 

In the following we will discuss these ideas in more details. In the first step we will show how simple linearizations 
will provide first ideas to the general problem. There we see how the nonlinearity changes the character of the dynamics 
completely, when the chains are strongly elongated. The simple analytic model will also provide the differences in 
the behavior of the two components R± (s) and i?| | (s) . In the second part of the paper we will use simulations to go 
beyond the linear analytical approach. We shall see that the simulation results provide important new information 
on the chain dynamics in the high stretching regime. 

II. GAUSSIAN SELF-CONSISTENT APPROACH 

To proceed in a first step with the analytical theory we have chosen to employ a self consistent Gaussian approach, 
i.e., a variational technique, which allows to get a first feeling for the behavior of such strongly elongated chains. 
The first step is, that for such stretched chains with finite extensibility we cannot argue anymore in looking at the 
isotropic chain variables themselves. Therefore the chain vectors need to be separated into components parallel and 
perpendicular to the stretching direction. Despite the mean field, and the linearization character, this will provide 
some important insight to the physics of the problem. 

The dynamics of a stretched polymer chain at high degree of stretching is essentially nonlinear. This is due to 
the necessity of large forces being applied over a short region of space so as to prevent the chain from stretching 
indefinitely. While in the Rouse (linear) model, the linear springs that provide the connectivity of the chain fail to 
meet this condition, it is well satisfied by the anharmonic springs, described by the frequently used finite extensible 
nonlinear elastic (FENE) potential. However, we treat first the relaxation dynamics of a stretched polymer chain by 
solving the nonlinear Langevin equation within the framework of the Gaussian self-consistent approach (GSC). To do 
so we considered a polymer chain with N monomers connected by FENE springs, with one chain end fixed at the origin 
of coordinates and pulled by force / at the opposite chain end. Due to the external force /, there emerges a preferred 
direction in space leading to symmetry breaking in the problem, caused by extending the polymer in direction of 
the acting force. This effect splits the dynamics of the chain into two directions^ with respect to monomers motion, 
R(s,t), one parallel to the force and the other perpendicular to it, leading to distinct longitudinal and transverse 
relaxation times. For the equation of motion of the statistically averaged position, Ru(s,t) = (-Rii(s, f)), in direction 
parallel to force, we have 

ggj|(M) / SV\R(s,t)-R(s-l,t)] \ / 6V[R(s,t)-R(s + l,t)] 
^° dt + \ SR\\(s,t) / + \ SR {] (s,t) 



\ - fSsN = (7) 



where V[R(s, t) — R(s ± 1, i)] refers to the non-linear potential (FENE in this case) which links the nearest neighbors 
along the chain backbone and £o denotes the friction coefficient of a single monomer. We introduce the notation, 
R»(s,t) = (R z (s,t)}. In direction perpendicular to force, the first moments are R x (s,t) — R y (s,t) = Rj_(s,t) = due 
to axial symmetry. Additionally, we present the equations of motion for time displaced correlation functions; in the 
perpendicular direction C±(s,n,t,t / ) = (R±(s,t)R±(n, t')) we write 

dC^.n,t,1>) + / R± t , SVms, t) - R(« - It)] + V[K(s, t) - R(, + l,t)]\ = {fltL( (g) 
ot \ oR±(s,t) I 

and for the parallel counterpart C\\(s,n,t,t') = (R\\(s,i)R\\(n,t')) 

+ (*. (», O ^M)-R( a -l,0] + mM)-R( a + M)] ^ _ ^ (B> ^ = (J?|| (n> 0fcj (S) t)) (9) 

where /i||(s,£) and /ij_(s, i) are random Gaussian (S-correlated forces in parallel and perpendicular directions respec- 
tively. We also study the equal time correlation functions C\\±(s,n,t,t) = C\\.±(s,n,t) which in the steady state 
regime do not depend on the time t. For the perpendicular direction, we have 

to dC ^ n,t) + K±(s, s - l)[C x (n, s, t)-C ± (n,s-l, t)} + K ± (s, s + l)[C±(n, s, t)-C±(n,a + l, t)} 
+K±(n,n - l)[C±{s,n,t) - C±(s,n - l,t)] + K±(n, n + l)[C x (s,n,t) - C±(s,n+l,t)] = 2S ns k B T (10) 

and for the parallel, 

dAu (s n t) 

Co — 11 „ ' + K\\ (s, s - 1) [An (n, s,t)-An(n,s-l,t)]+Kn(s,s + 1) [A n (n, s,t)-A n (n,s + l, t)] 

ot " " " " " " 

+K\\(n,n - l)[A\\(s,n,t) - A\\(s, n - l,t)] + K\\(n, n + l)[A||(s,n, t) - A\\(s,n+ l,t)] = 25 ns k B T (11) 

Here, we define A\ \ (s, n, t) = C\\(s,n,t) — Rn(s,t)Ru(n, t), and the effective nearest neighbor constant, K\\ j_(s,s± 1), 
is defined by eq. (|A3j) in the Appendix. The relationship (2Rx(n, t)h±(s,t)) — 2<5„ s fcsT is also derived explicitly in 
Appendix O 



A. First cumulant expansion 



Consider first the term ^ £KI5fe*2. — b*)l\ which is the average force experienced by each monomer. The calcu- 
lation yields 

'5V[R(s,t) - R(s- l,t)] x 



SR\\(s,t) 



6R\\(s,t) 



d 6 rS[H(s, t) - R(s -l,t)- r]V(r 



\ dr^ R ( s ' ^ ~ R ^ s ~ lf *) ~ r ] ) V ^ = / d r (W s > *) - R ( s - 1. *) - r ]> 



ay(r) 
dm 



The key point is then to determine the difference (<5[R(s, t) — R(s — 1, t) — r}). This is accomplished by employing 
the following crucial approximation (<5[R(s,t) — R(s — l,t) — r}) sa J[(R(s,i) — R(s — l,t) — r)] which comes from 
retaining only the first term in the cumulant expansion. Taking into account that V(f) = — |fci?&Qln(l — p-) and 



= r» + r\ one arrives at 



/ 



d 3 r(S[R{s,t) -R(a- l,t) - r]) 



9y(r) , fl||(a,t)-B||(«-l,t) 



6o 



B. Steady state solution 

In order to obtain the steady state solution of eq.0, we set dR ^' t * > = 0. Then, we can find -Rii(l) by adding up 
the N terms for the -Ry(s), requiring that the boundary conditions of the problem are Ru (0) = (first monomer is 
fixed at the origin of coordinates), and R\\(N + 1) — Rii (N) — (free end condition). 



Eventually we get 



from which one can easily derive i?||(l) 
function of /. 



kpb\ 
~2j 



(12) 



k ^— 1 + «/l + -gf^j- In Figure JT]) we plot the elongation i?||(l) as a 
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FIG. 1: Variation of the bond length between adjacent beads in a chain with applied tensile force /. Eventually, for large forces 
the bond length approaches asymptotically its maximum value bo- 

Applying an iteration procedure, a general expression for R\\(s) can be found, R\\(s) = S-Rii(l), 1 < s < N. 

Setting dCh g/'" = 0, the steady state solution for the equal time correlators can also be obtained. Then, for the 
perpendicular direction, the final equation reads 

K_l(s, s - l)[Cx(n, s) - C±(n, a - 1)] + K ± (s, s + l)[C±{n, s) - C±(n, s + l)] 
+K ± (n,n- l)[C±(s, n) - C_l(s, n-l)] + K ± (n,n + l)[C ± (s, n) - C±(a,n + 1)] = 25 ns k B T, (13) 

and for the parallel 

K\\ (s, s - l)[A\\ (n, s) - A\\ (n, s - 1)] + K\\ (s,s+ 1)[A {1 (n, s, t) - A\\ (n, s + l)] 
+K ll (n,n-l)[A {] (s,n) - Ay(s,n- 1)] + K\\(n,n + l)L4||(s,ra) - A ll (s,n+ 1)] = 25 ns k B T (14) 

The solution for the time- independent correlator Cj_(s, n) is found by rearranging eq. (|13p which makes it possible 
to take its continuum limit and rewrite it in a more convenient way. Eventually, we get 



K i 



<9 2 Ci(n,s) a 2 C ± (n,s) 



= -2rJ(n - s)fc s T 



(15) 



ds 2 dn 2 

where K±(s, s± 1) = JCx(n, n± 1) = K± = — / ^ F (i) . ^ . It is better here to define this constant by the relative degree 



f>0 



of stretching A = and its maximum value, X max = = V~N (where Nbo denotes the total contour length 

of the chain) . Thus one can rewrite K± = — -r-^ — ^ • 

Equation (fl~5|) is a linear equation which can be easily solved for example by means of the Green's function method^. 
In the same way, one can rearrange equation (|14[) and proceed in a similar fashion. The final solutions reads 



4 k B T 



2N 2 



^[NK l± (2p-l)H 2 



(2p-l)im\ /(2p-1)tts 
2N ) Sm { 2N 



(16) 



where K\\(s, s ± 1) = K\\(n,n ± 1) = K\\ = kpj — V A„ lax j ^ ^ These two coupling (effective spring) constants behave 

therefore differently in the two directions and depend in a detailed way on the deformation ratio. Both spring constants 
show a singularity at the maximum deformation. We will use these results again in the last section of these paper, 
when we discuss the limiting cases and the physics of the numerical results, which will follow later. 



C. Time displaced correlators 

The steady state solution for C\ix(n, s) offers an initial condition from which the time dependent problem can be 
solved. Following the GSC method, we are now in a position to solve the equations of motion also for C±(s,n,t,t'). 
Using the causality conditions of physical systems (R± (n, t')h± (s, t)) — for t' < t, and the same calculations outlined 
in the previous section, after taking the continuum limit, one finally arrives at 

dC ± (n,s,t,t') f d 2 C ± (n,s,t,t') \ 

^ — di — =M — a? — )■ (17) 

This is one-dimensional diffusion equation that can easily be solved to give 



(18) 



where f3 p = < ' 2p ^\ ~ and we have chosen t' = without loss of generality. The amplitudes A p may be determined 
from the initial condition since C±(n, s, 0, 0) = C±(n, s). Equating the expression for C±(n, s, t, 0) at t = 0, (eq. fl"8|) . 
and Cj_(n, s), one arrives at 

4 k B T 2N 2 (2p - l)im , , 

Ap =N-Kl (2 P -l)^ Sm 2N ■ (19) 

Evidently, A p is a function of (p, n) since C±(n, s, t, 0) is a function of (n, s, t). Now, it is possible to write down the 
final expression for C±(n, s, t, 0) as 

4 k B T S pm ^-p B *jL t ( (2p - l)im\ ( (2m - l)irs 

N K, T2pl) 2 ^ 2 i (2m^lj 
D =l m=l 2JV 2 1 



Cj_(n, s,t,Q) = > > —-7^-77^ — 2 w TT2-^" e 5(1 sm sin — (20) 



Following the same reasoning one can obtain an identical solution for the parallel direction too by replacing A"j_ — > Ku . 
From the last equation, the (Rouse) mode amplitudes can be calculated in a straightforward manner: 

(X pU{ ( t)XmUl (0)> = M.^^— exp (W^ft) = ^S Pm ex P (-^) , (21) 

where T p ±^ = j^' 11 = (2 P -\yiJ^K ± ' * s ^ ne perpendicular (parallel) relaxation time and kpx,\\ = ^ P 27V ~ anc ^ 

One should stress that the relaxation time decreases in both directions as one increases the degree of stretching 
since r oc 1 — i^j^- — ) • This behavior can be understood, if we take into account that the stretching of the chain 

makes the springs effectively stiffer. 

We must also point out that within the present approach the chain behaves effectively as Gaussian, r oc N 2 p~ 2 , 
which is plausible since chain stretching rapidly reduces the role of the excluded volume interactions between the 
beads. This feature reveals the nature of the approximation that we have used (see above) which essentially drives 
out the coupling between different modes but still retains the physical intuition of a zero relaxation time for maximum 
stretching. 



III. SIMULATION RESULTS 



To test the validity of the approximations made in the previous section, we performed computer simulations. Two 
different numerical schemes have been chosen to test the relaxation dynamics of stretched polymer chains and also to 



test the validity of analytic predictions: Monte-Carlo and Molecular Dynamics. Physical intuition suggests that for 
highly stretched polymer chains the two methods are expected to give different results. This difference stems from 
whether the inertia term in the simulations is neglected or not, which for large extensions of the polymer plays a 
decisive role. Since Monte-Carlo simulations do not consider mass in the model, the corresponding molecular motion 
of the chain must be of an overdamped nature (oscillations are not possible). In contrast, Molecular Dynamics takes 
inertia into account and it is possible to pass from an overdamped motion (small extensions) to an underdamped one 
(large extensions), in which oscillations can occur. In what follows we analyze the simulation results derived by both 
methods, and compare them to the predictions of the previous section. 



A. Monte-Carlo simulations 



We perform the simulations with an off-lattice coarse-grained bead-spring model which has been frequently used 
before for simulation of polymers^. Therefore here we will only describe briefly its relevant features. The effective 
bonded interaction between nears-neighbor monomers is described by the FENE (finitely extensible nonlinear elastic) 
potential. 



Ufene = -K(l - l ) 2 ln 



I -la 



(22) 



with elastic constant K = 20. The maximum bond extension is l ma x = 1, the mean bond length Iq = 0.7, and the 
closest distance between neighbors l m in = 0.4 

The nonbonded interactions between monomers are described by the Morse potential. 

Um ^ T = exp(-2a(r - r mm )) - 2exp(-a(r - r mm )) (23) 

with a = 24, r mm = 0.8, eM/ksT — 1.0. 

The size of the box is 64 x 64 x 64 and the chain length N = 32. Of course, a larger number of monomers could in 
principle be used but this is not necessary since the most relevant features of the solution can be well illustrated with 
this number of monomers. The standard Metropolis algorithm was employed to govern the moves with self avoidance 
automatically incorporated in the potentials. In each Monte Carlo update, a monomer is chosen at random and a 
random displacement is attempted with Ax ,Ay ,Az chosen uniformly from the interval —0.5 < Ax, Ay, Az < 
0.5. The transition probability for the attempted move is calculated from the change AU of the potential energy 
U = Ufene + Um 'as W = exp(— AU/ksT). As usual for the standard Metropolis algorithm, the attempted move 
is accepted if W exceeds a random number uniformly distributed in the interval [0, 1). We employ 10 4 runs of 2 16 
Monte-Carlo steps in each program run. 



1. MC Simulation results 



During the simulation we compute ry and t_l, and the time displaced correlation functions between different modes, 
C pq (t) = (X p n ±(t)X q i\ j_(0)) (where p, q label the mode numbers) as a function of the relative degree of stretching 
X/Xmax- Here, the average is taken over different intervals of time t. 

It is straightforward to compute also theoretically this quantity in the case of Gaussian chains without excluded 
volume interactions^. This leads us to 



C pq (t) = C8 pq exp — ) p= I..N 



-t 



(24) 



where C is its initial value. For the relaxation time t p we obtain 



N 2 £ b 



3fc B T(p- 1/2) 2 tt 2 



(25) 



which clearly shows the t p cx N 2 p~ 2 dependence of the relaxation time due to the harmonic nature of the chain 
interactions. Obviously, there is no p = translational mode since the chain, fixed at both ends, cannot diffuse. For 
the highly stretched case it is not clear how the chain will behave, since the finite extensibility plays a fundamental 
role and the forces along the chain are non- linear. 



We perform simulations for stretching degrees of 0.5 < X/X ma x < 1 which are shown in Figure [2J One can 
see an almost perfect agreement between the theoretically predicted (GSC) and simulated (Monte-Carlo) values of 
the relaxation time riy of the first mode in parallel direction. For the perpendicular direction the matching is not 
so good but there still exists a qualitative agreement between analytical and numerical results (the discrepancy 
at X/Xmax — > 1 is due to numerics since the acceptance rate of the elementary displacements depends strongly of 
the attempted jump distance Ax, Ay, Az). The relevant feature here is that the relaxation time in both directions 



5000 rr 




FIG. 2: Monte-Carlo simulations of the first mode relaxation time for the perpendicular rj_ and parallel direction tii (MC 
steps) as a function of relative degree of stretching 0.5 < \/\ m ax < 1. Comparison with analytical (AN) and numerical results 
(MC). 

decreases as the chain is gradually stretched to its final limit. Evidently, the more stretched the chain is, the "stiffcr" 
it becomes and, consequently, the relaxation time goes down. There is also a pronounced difference between the 
parallel and perpendicular relaxation times for the same X/X max . The correlations in parallel direction decay faster 
(shorter relaxation time) than in perpendicular direction. However, this difference tends to disappear as the degree 
of stretching increased. 

An important question concerns the degree to which the stretching affects the correlation between the first modes. 
In Figure [3] (a) and (b) we plot the autocorrelation (-X^m j_(0) 2 ) and cross correlation (^i|L_L(0)X m ii j_(0)) functions 
as regards the mode number m for different degrees of X/X ma x- We want to emphasize two features. First, one should 
note that the power law of r p ii j_ oc p -1 - 98 is observed for all extensions and in both directions as this can be clearly 
seen from the insets of Fig. IS (a) and (b). This behavior corresponds to our theoretical predictions and reflects the 
fact that the chain behaves like Gaussian at least as long as 0.5 < X/X max < 0.975. 

Second, we find that the modes of the same direction seem to be orthogonal to each other for 1 < p < 31, as in the 
linear (Rouse) case. For the cross correlation terms between different directions, (Xi\\(0)X p ±(Q)) , this is not the case 
at least for the lower modes 1 < p < 10 which suggests that finite extensibility couples these modes in a certain way, 
see Fig. |U 

In conclusion, one may claim that this coupling effect, due to the nonlinearity of the interactions, seems to be stronger 
between modes of different directions and has apparently no effect on modes in the same direction. Generally, only 
the first modes (1 < p < 10) appear to couple since for modes with p > 10 (that is, for short distance motions) the 
chain does not "feel" the fixed ends and the behavior is that of Rouse dynamics. 

B. Molecular Dynamics results 

Although real polymers do have mass, the inertia term is almost irrelevant as compared with the friction term, at 
least as long as the chain is not extended,. For high extensions, however, due to the large accelerations that each 
monomer is subject to, the inertia term is comparable and even larger than the friction term and the whole chain 
behaves as a string under tension in which oscillations are no more overdamped. 

In this section we check the relevance of the inertia term, which should be important at the high stretching limit, 
and examine its physical consequences. To this end we run Molecular Dynamics simulations in which the mass of 
each monomer is accounted for in the corresponding equations of motion. For the implementation of the method we 
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FIG. 3: (a) (0)X p || (0)} against mode number p for different relative extensions A/A max . In the inset we plot the p 
dependence of the function (X p ^ (0) 2 }, we show that the slope is almost equal to the predicted theoretical value oc p~ 2 (b) The 
same as in (a) for the perpendicular direction. 
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FIG. 4: (X^ (0)Xpx(0)} as a function of the mode number p for different extensions X/X r . 



employ a successful numerical scheme developed and proved earlier by Dimitrov et al.— . The model is again that of a 
simple coarse-grained bead-spring chain, originally proposed by Kremer and Grest- 19 , which has been widely and very 
successfully used for MD simulations of polymers in various context o 20 ! 21 . Effective monomers along the chain are 
bound together by a combination of the FENE potential, Eq. (j2"2"j) . where K — 15 and the maximum bond extension is 
Imax = 1.5(7, and a Lennard- Jones potential, responsible for the excluded volume interaction between the monomers. 
a = 1 is the range parameter of this purely repulsive Lennard- Jones (LJ) potential which is truncated and shifted to 
zero in its minimum and acts between any pairs of monomers. 



U LJ (r) = Ae LJ [(a/r) 12 - {a/rf + l] ,r < r c = 2 1 ' & a 



(26) 



The parameter e^j, characterizing the strength of this potential, is chosen unity. Molecular Dynamics (MD) simula- 
tions were performed using the standard Velocity- Verlet algorithm 22 , performing typically 1.5 x 10 9 time steps with 
an integration time step St = O.Olio where the MD time unit (t. u.) to = (rj 2 m/48eLj) 1 / 2 — 1/V48, choosing the 
monomer mass = 1. The temperature was held constant by means of a standard Langevin thermostat with a 
friction constant £o = 0.5 



1. MD Simulation results 



To make a comparison between Monte-Carlo and Molecular Dynamics simulations, we study the same cases with 
the same polymer chain of length N = 32. This time, 2 s MD steps in each program running have been performed. 

Typical data for the time, equal- mode correlation functions C pp (t) — (X p \i ±(t)X p \^±(0)} are presented for a partic- 
ular value of X/Xmax — 0.83 and for the mode p — 3. Figures [HK andjSb confirm the importance of taking inertia into 
account, manifested by the oscillations observed in the correlation functions. In Fig. [SK, which corresponds to the 
mode in parallel direction, one can recognize at least three frequencies in the correlation function, exposed by smaller 
satellite peaks in the Fourier spectrum in the inset. This broad-band spectrum does not mean that only 3 frequencies 
are present in the correlation function. One may rather claim that a continuum of frequencies of decreasing amplitude 
appears instead, reflecting the fact that these frequencies are clearly discernable from their neighbor values. Appar- 
ently, this effect manifests itself less visibly for the perpendicular direction, at least for this value of X/X max where 
only one peak can be well distinguished (see Fig. [SJs). This single frequency appearing in the correlation function 
allows us to consider it as a linear underdamped motion of the chain for this case. 




FIG. 5: (a) (X p u (t)X p ±(0)) versus time (in MD time units) for p = 3. In the inset the Fourier spectrum of the same function 
shows a sharp maximum, corresponding to the high basic frequency. Two smaller satellite peaks correspond to lower frequencies, 
also present in the mode dynamics, (b) The same as in (a) but in perpendicular direction. The Fourier spectrum reveals the 
presence of a single well-defined frequency in contrast of other two, less visible, smaller peaks. 

At this stage we can ask for the physical meaning of these satellite peaks shown in the above figure. To answer 
this question, we plot together, in Figures [6^ and [6b the Fourier spectra (frequencies) of several mode correlation 
functions at three different degrees of extension. The first five modes in parallel direction are shown in Fig. [5^ 
for X/X max = 0.512 ,0.704 ,0.953. One can readily see that with growing stretching, X/X max , a central peak, 
corresponding to a principal frequency of oscillation, is formed. However, this is not the case for X/X max = 0.512, 
where such a peak is absent and the mode dynamics is spread almost uniformly over frequencies lower that 10~ 4 . 
In the other two cases X/X max = 0.704 ,0.953 this spreading disappears and narrower peaked functions come out, 
suggesting that the mode motion is realized essentially by means of a single frequency. Now, we claim attention to 
a particular feature of Figured^. Take for example the case X/X max — 0.704 for p = 1. There, by a vertical line 
we demonstrate that the frequency corresponding to this mode (p = 1) coincides exactly with the satellite frequency 
peaks which appear in the higher mode correlation functions. Looking deeper, we can see that this is not prerogative 
of the first mode and is observed also for p > 1. It takes place for larger extensions too. Clearly, this shows that the 
satellite peaks are exactly the main frequencies of higher mode correlation functions. At the same time it indicates 
coupling of lower modes with higher ones which is accompanied by exchange of energy between them. Evidently, 
this coupling effect becomes more significant with growing mode number and, most notably, at lager extensions (see 
curves for X/X max — 0.953). 

Figure [SJd shows the same quantities for the modes in perpendicular direction. In this case the situation is similar 
whereby the satellite peaks indicate that coupling between the modes is effectively weaker, and may be observed at 



the largest extensions only. 




FIG. 6: (a) Fourier spectra of correlation functions {X p ^(t)X p ^(0)) for the first 5 modes in parallel direction at three degrees of 
relative extension X/X ma x- The dash-dotted line is plotted as a guide to the eye. (b) The same as in (a) for the perpendicular 
direction. 




FIG. 7: (a) Idem Figure [S] where a linear scale was used to show the broadening of the half-width of the main peak in the 
Fourier spectra of correlation functions (X p u(t)X p «(0)) for large values of X/X m ax = 0.704 0.953. (b) The same as in (a) but 
this time the broadening if exists, does not appear to be significant. 

More detailed information of the mode coupling effect can be extracted if one plots, for higher extensions X/X max — 
0.704 , 0.953, the same functions as in Figured This time, a normal scale in the frequency axes (abscissas) is employed 
to facilitate visualization. In Figure [7^,, two significant effects can be well distinguished. On the one hand, we can 
appreciate a slight broadening of the half-width of the main peak with increasing extension for the same mode. On 
the other hand, a similar broadening occurs with increasing mode number, but this time at constant extension. Also 
the height of the main peak decreases with growing p. From a physical viewpoint, a broad-band spectra centered 
in a given frequency is associated with a continuum of frequencies (as mentioned above) around this single value. 
Returning to our problem, this means that as soon as we increase extension, more and more frequencies appear to 
make larger contributions in the dynamics of single mode correlation functions. A distinct situation accounts for the 
perpendicular direction. Noticeably, this broadening is not seen in Figur^Zb, and the peak height remains constant 
as one may verify from the similar graph depicted. 



It is interesting to analyze the behavior of the principal frequency lu, corresponding to the peak in the Fourier 
spectrum regarding the polymer chain extension. 

In Figure EH we demonstrate that the characteristic frequency in parallel direction u\\ grows with increasing ex- 
tension. At fixed value of the elongation X/X max , Wm systematically increases with growing mode number p. For a 
special value of X/X max rs 0.64, however, one observes a local maximum in which becomes increasingly pronounced 
as the mode number p increases. In the inset to Figure EH, where we show the variation of the relaxation times 
t p of the modes with changing elongation X/X max , the region around X/X max 0.64 marks the end of the steady 
decrease of t p with growing degree of stretching. Thus, it appears that at some threshold extension X/X max w 0.64 
the chain dynamics undergoes a qualitative transition from an overdamped motion, determined essentially by friction, 
to ballistic motion governed by inertial effects. 

In the overdamped regime the MD results are qualitatively consistent with the analytical predictions and MC data, 
indicating a clear decay of the relaxation times with increased stretching. Beyond some critical stretching, however, 
the dynamics changes and the polymer chain behaves increasingly like a string under tension. Then oscillations rather 
than random displacements of the monomers become important and the relaxation times of the different modes start 
to increase once again over a short range of chain elongations. Upon further stretching these relaxation times stay 
nearly constant. Eventually, in a regime of very strong stretching for X/X max > 0.8 the system becomes strongly 
nonlinear in nature and it becomes difficult to determine a well defined single "relaxation time" . 
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FIG. 8: (a) Fourier spectrum of parallel mode correlation functions, (X p \\ (t)X q n (0)) against chain extension. In the inset one 
shows the variation of relaxation times t p for the first five modes with chain elongation X/X ma x- (b) The same as in (a) for the 
perpendicular direction. 

Figure shows similar plots for the perpendicular direction. Rather than a maximum, the threshold position at 
X/X max ~ 0.64 marks the onset of a stepper increase of the main frequency coj_ with growing elongation whereby 
at X/X max — > 1 this frequency appears to diverge lu±_ — > oo! The relaxation times tj_ of the various modes in 
perpendicular direction behave differently as compared to their counterparts of Fig. EH. While in the overdamped 
regime at X/X max < 0.64 they decay steadily with growing elongation, in the inertial regime here after the threshold 
they merge to a nearly single value which remains largely unchanged for X/X ma x > 0.64. Remarkably, since the 
oscillations of the string in perpendicular direction persist up to very large extensions, X/X max ~ I, a well defined 
relaxation time can be found. 



2. Approximate analytical model 



It is possible to explain the observed behavior of the principal frequencies in the Fourier spectra of the mode-mode 
correlation functions as well as that of the relaxation times if one considers the linearized equations of motion |7|. 
The intention here is to capture the essential features of the main frequencies wim dependence on mode number p 
and understand the insensitivity of the relaxation times with regard to chain elongation. 

Our model of a chain pulled by force / and fixed at the origin allow us to think of the chain as being fixed in space 



with elongation R»(N). One can then write an approximate equation of motion for R(s, t) using the information 
provided by the GSC method of section HT1 

First, we will decompose the position vector as R(s, t) — R(s, t) + (R(s, t)) 1 where R(s, t) is a new vector measuring 
the deviations of the position vector from its average position (R(s,i)). Then, making use of previously defined 
quantities and expression (1), we can propose the following linearized equations 

d 2 R j {s,t)\ /dRj(s,t)\ /d 2 R 3 (s,t) 



j£r ) + ^\dr L /~ Kj \ ds 2 - )-f S ^ = ° 3=x, y, z (27) 

where we have replaced the exact force by a linearized (approximate) elastic force with elastic constants Kj = Ku , K±, 
obtained in section Hi Bl 



/ SV[K(s,t)-K(s-l,t)} \ / SV[R(s,t)-K( s + l,t)} \] ^ I d*R lu± (s,t) \ 

Separating the problem in the parallel and perpendicular direction as before and using the previously defined decom- 
position of vector R(s, t), since equations (|27|) are linear, we finally arrive at 

mb \ at* I + e ° \ at / ~ K ^ \ b / = (29) 

Introducing the Fourier expansion for (R^ j_(s,t)) = X^>(^p||,-i-W) sm (*7jr) ano - inserting it in (|29|) . we obtain after 
multiplying it by ,j_ (0) 

J2 (m b C plU ±(t) + ZoC plu (t) + K\ U± (^) 2 C plu (t)\ sin (^) = (30) 

p \ / 

where C p \\ ±(t) = dC p \\±(t)/dt and C p u±(i) = j_(t)X p ii j_(G)). To satisfy f[30|). we require that 

m b d A ,x + e C p || , ± + Aj, , ± (^) 2 C P || , ± = (31) 
Assuming that C7p|| ,x W = Cpii x ex P ( — Tp||,_L^)j we finally obtain the eigenvalues of the linearized problem 



lpl± ~ 2m b ^ 4m 2 m b \n) [6 ' 

Clearly, these eigenvalues are real, showing a purely exponential decay, or complex, indicating damped oscillations, 
depending on the value of the radical in (f3"2"|) . Hence, these cases must be treated separately. 

• Overdamped case 

In this case, the radical is real, provided 

(33) 

2m b J m b \NJ k ' 

and it is possible to identify 7 p ||,x with the inverse of the relaxation time Tpii x = V r p||,-L- 

Returning to (|32[) and retaining only the first term in the Taylor expansion of the radical, vl — x 2 sa 1 — ^£ 2 , we 
find two values of 7pii x or T p||,J- : 



pII.-L pII.J- 



(34) 



since eq. (|3"Tj) is of second order in time. Evidently, to recover the theoretical result of section [TT] for the relaxation 
time Tpi^xj one must neglect the inertia term (inserting, for example, m b = 0) in (|27[) whereby only t^x survives. 



• Under damped case 



In this case the radical is negative, 




Now, 7 p ii j_'s are complex quantities meaning that C p \^j_(t) corresponds to damped oscillatory motion. Clearly, 
the real part of 7 p ||,j_ describes the degree of damping and the imaginary one gives the oscillation frequency which 
can be called If we still consider the real part of 7 p ||.j_ as the inverse of a new relaxation time r p » j_ for the 

underdamped case, we find that Tpii j_ is equal to 




FIG. 9: (a) Main frequency of the parallel mode correlation function u> p « as a function of mode number p for different extensions 
in the underdamped regime. In the legend X/Xmax means analytical approximation, MD , Molecular dynamics simulations, (b) 
The same as in (a) in the perpendicular direction. 



where i = is the imaginary unit. This can be approximated for large p or large extensions 3> 1) by 

1 ^ and c p|u « ± f^V /2 m (37) 



7jj| iX 2m b \ m b J \NJ 

where a Taylor expansion of the radical up to first order was carried out. 

The last equation reveals some remarkable features. First of all, we can see a linear dependence on p of the main 
frequencies in the mode correlation functions w p |ij_. Secondly, an increasing of the magnitude of these frequencies 
upon increased stretching is also demonstrated. Figures^ and b support these observations showing the reliability 
of the approximations. We plot and uj p ± [SJs against mode number p for different extensions. Note that these 

findings apply only to the underdamped regime, X/X max > 0.64. 

Moreover, a nearly constant relaxation time, T p u ± , which doesn't depend on mode number or extension, comes out 
from equation ([57]) too. Upon qualitative comparison between Figures [5£i and b (inset plots), we can conclude that 
for the perpendicular direction [SJd, where a single relaxation time can be well defined even at very high extensions 
X/X max 0.95 the agreement is good. Instead, for the parallel direction [SK, a constant value of the relaxation time 
agrees with simulation results up to the extension where a single relaxation time can be well defined. 



IV. CONCLUSIONS 



In the present paper we have studied the dynamics of a self-avoiding polymer chain at different degrees of stretching. 
Analytical work as well as computer simulations provide a rather consistent picture of the qualitative changes which 



the polymer dynamics undergoes upon gradual increase of the degree of stretching. In general, one observes two 
regimes of chain dynamics, depending on the degree of chain extension. In the first one, that of the friction dominated 
overdamped motion of the monomers, both analytic predictions as well as Monte Carlo results suggest a consistent 
picture of relaxation time decrease with growing stretching of the chain up to a threshold value X/X ma x ~ 0.64 
whereby the relaxation time parallel to stretching, T|| is always considerably smaller (about one half) of that in 
direction perpendicular to stretching, t±. For T|| the agreement between analytic results from the GSC approximation 
and MC data is perfect on a quantitative level whereas for t± it is at most qualitative. As expected, the MC results 
for the relaxation time versus mode index relationship yield t p cx p~ 2 . 

A transition regime at \/\ ma x ~ 0.64 where friction and inertia terms are of the same order of magnitude separates 
the overdamped regime from an incrtial regime at higher degrees of stretching. This latter regime is strongly nonlinear 
in nature and a faithful description of the polymer dynamics can be produced by means of a MD simulation as well 
as by an approximate linearized analytical model. 

Among the most salient features of polymer dynamics in this second regime of strong stretching, we find that 
normal modes are coupled to each other and can interchange energy, as expected in a system with strongly anharmonic 
interactions. These anharmonic contributions to the bond (FENE)-potential come into play at strong stretching only 
whereas in a weakly extended chain the monomers make the bonds oscillate around the equilibrium bond length by 
means of nearly harmonic forces and the modes are largely independent. 

The Fourier analysis reveals the existence of principal characteristic frequencies of oscillatory motion for each normal 
mode. Along with this feature, a broad-band spectrum suggest the existence of a continuum of frequencies around the 
single peaks, which apparently make a non-vanishing contribution to the total motion. These principal characteristic 
frequencies have been shown by an approximate linear analytical model to scale linearly with mode index p. The 
frequencies also grow upon extension, as expected. A slight broadening of the half-width of the main frequency peak 
with increasing extension as well as with increasing mode number at constant extension, is revealed for the parallel 
direction. This means that as soon as we increase extension, more and more frequencies appear to contribute in the 
polymer dynamics, reflected by the single mode correlation functions. This broadening is not evidenced, at least up 
to X/Xmax = 0.95, for the perpendicular direction. 

The relaxation times of the modes depend significantly on the considered direction. While in the parallel direction 
the relaxation times increase strongly in the transition region but then the possibility to define a single relaxation 
time for too large extensions X/X max > 0.8 becomes questionable, in the perpendicular direction the relaxation times 
remain more or less constants up to X/X max « 0.95. It is remarkable that these results are also confirmed qualitatively 
by our linearized approach. Of course, one must keep in mind that such large extensions of a polymer chain might be 
of academic interest since bond rupture make take place if one approaches the tensile strength of the macromolecule. 
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APPENDIX A 



In calculating terms like ( ^^[R-tM) Mj] \ one neec i s to tackle the term (<5[R(s, t) — R(s — 1, t) — r]) - see Section 



III Al To this end we make the following approximation 

(S[R(s, t) - R(s - 1, t) - r]) w 5[(R(s, i) - R(s - 1, i) - r> 
Taking into account that V{r) = — ^kph^hxil — pr) and r 2 = rj + r\ 1 one can set 



(Al) 



d 3 r(S[R,{s,t) -R(s- l,t) -r])-^ » fc F 11 V ; 11 V 



Hy(a,t)-JJ|,(8-l,t) 



which is the final expression that must be inserted in the equation of motion. 
The calculation of Ay (s, n, t) = Cm (s, n, t) — Ru (s, t)Ru (n, t), proceeds as follows: 

From equation (J9j) we have to calculate terms like ( Ru (n, t) IZIBii^z^^i^i*)] \ Using the 5-function, one can write 



R\\{n,t 



SV[K(s,t)-IL(s-l,t)] 

SR\\{s,t) 



R\\(n,t) / d 3 r(S[R(s,t) - R(s - l,i) - r]> 



(A2) 



Applying Wick's theorem 2 ^ we have 

JV[R(s,t) - R(s - 1, i) 



R\\{n,t)- 



)R ( f) ^ = (R l{ (n,t)) J d 3 r(8[R(s,t) -K(s - l,t) - r])^^- + ((R }] (n,t)R }] (s,t)) 

-iR^n,t))(R l{ (s,t))) ( — S - ^ J d 3 rS[R(s,t) - R(« - - r]^^ + ((-Ry (n, t)fl y (a - l,f)> 

,^,( S -M)/ d3r5[R(M) ' R(s - M) " r] ^r) 

<9 2 V(r) 



>( 57?||(s,i) 
-<«„ (f», *)> <-R„ (« - 1, t)» 



7?H (n, t) J d 3 r(5[R(s, t) - R(s - 1, i) - r]) ^^ + [A,, (n, s, t) - A\\ (n, s - 1, t)] ^ d 3 r(<5[R(s, i) - R(s - 1, t) - r]) 



drf. 



where the rhs was obtained after integrating by parts. For the term (^Rj_(n, ty ^^jR^j^ 1 ' t ^ ^ ! one h as to replace 

the parallel variables by the perpendicular counterpart, taking into account that R±(s) = 0. 
Concerning the meaning of Kuj_(s 1 s ± 1), we define an effective elastic constant: 



K lA _(s,s±l) = J d 3 r(5(R(s,t) -R(s±l,f) - r)) 



d 2 V"(r) 
9r 2 i 



(A3) 



Applying the same approximation as before, and after taking the second derivatives of the potential V(r) with respect 
to perpendicular and parallel directions, we have for the perpendicular effective constant 



K ± {s,s±l) = k F - 



1 



and for the parallel one 



K\\{s,s±l) = k F 



1 - 



\R\\{s,t) 


--R||(s±l,t) 


1 2 




bo 




'Ry(s.t) 


-R n (s±l,t)~ 


2 


bo 




\R\\{s,t)- 


-H||(s±l,i)l 


2\ 2 







(A4) 



(A5) 



The expression for the term (R±(n,t)h±(s,t)) can also be calculated by making use of the generalized Wick's 
theorem 2 ^: 



/ 5Rii±(n, t) 



(A6) 

As can be seen in—, the term ^ Jhjf^J^ ^ = 2^$ ns so tnat tne previous term gives 

{R\ ]t± (n,t)h\ U ±{s,t)) =5 ns k B T (A7) 
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